ICFEP/ABAQUS-UQ
Consisting of:
Enriched experimental of design + Bayesian inference
clearvars; clc; close all;
Part 1: Setup
Directory setup
- Set up working directory/input_file directory/output_file directory/subroutine directory
- Set up UQLAB core directory
- Set up working path (working path can be only under the file RUN)
%% Add UQlab core and subroutine directory
addpath("E:\OneDrive - Imperial College London\PhD_Work\UQlab\UQLab_Rel2.1.0\UQLab_Rel2.1.0\core");%UQLAB core
addpath("E:\OneDrive - Imperial College London\PhD_Work\SBI\Main\subroutines");%subroutine
% tell whether in the running folder
if ~contains(pwd(),caseInsensitivePattern('run'))
error('Working path is wrong!!!');
end
% define a working directory structrure for input/output
% decide whether reserve current FE results: true <-- delete; false <-- reserve
directory = Def_Dirs('true');
%% Initialize UQlab (not printing in the command window)
seed = 4;
rng(seed,'twister');
evalc('uqlab;');
Probabilistic input for EoD/parameters/discrepancy setup
%% Read input EoD parameters and observation error discrepancy hyperparameters range from csv
[Prior_forward,Prior_hyperdiscrep,Prior_EoD,AnParam] = readPrior(directory);
Observation data setup
% if this is a systhetic problem, add Gaussian noise percentage
Noise.scalar = 0;% Value as a percentage (e.g. 0.01-1.0)
%% read observation (ground truth) data
[observationcoord,observationsarray,gtarray,Noise,observation_stageTime,AnParam] = readobservation(directory,Noise,AnParam);
Analysis parameters setup
- Specify the number of seed/output fields/FE runs/Number of stage /enriched stage/parameters/percentage of training data over whole data
- Create "AnParam" structure to store the analysis parameters.
AnParam.FE_package = 'ICFEP';% choose FE package
AnParam.EoDseed = seed; %restore the seed
AnParam.N_RUN_EoD = 30; %number of FE runs required for EoD
AnParam.N_FE_stage = 70; %number of FE loading stages
AnParam.N_enrich = 71; %hope to enrich N_enrich-th stage for surrogate construction
% (1≤ AnParam.N_enrich ≤
% AnParam.N_FE_stage, if not in this range, no enrich)
AnParam.Switch_Predictive = "on"; % switch for prior/posterior predictive for real outputs
AnParam.TestDataRun = 5; %Test FE runs number
AnParam.TrainDataPerc = 1; %rangePerc = (e.g., 0.8-1), TrainRun = rangePerc * (Allrun - TestDataRun);
AnParam.Switch_surrogatePlot = "on"; % switch to show the test error/reconstruction error (only show in mlx)
AnParam.Switch_corner_scatter = "on"; % switch to show the posterior corner figures
AnParam.Switch_Sobol = "on"; % switch to show the posterior corner figures
AnParam.Switch_threeSigma = "on"; % switch to show the posterior paramters 3 sigma with stages
AnParam.Switch_RMSE = "1D"; % switch to show the posterior paramters mean and std with stages
AnParam.Trim = 1:size(observationcoord{1},2); %Trim range for FE output after interpolation
AnParam.DR = "on"; % choices for surrogates
AnParam.Surrogate = "PCE"; % choices for surrogates
AnParam.Switch_priorBound = "on"; % switch for prior forward boundary
AnParam.inferredCopula = "off"; % this button always encounter error with AnParam.Switch_priorBound
AnParam.figureDisplay = "Show-in-mlx"; % switch for display/save figs
AnParam.LL_Normalize = "off"; % normalize LL by MAP values
AnParam.LL_MAP_test = "off"; % test LL evaluation values at MAP and save in a csv file
% Used in readoutput.m Very difficult to get rid of this as .m is hidden by UQlab.
if strcmp(AnParam.FE_package,'ICFEP')
%%% Define global variables
global FE_ID
FE_ID = 0;
end
UQlink setup
%% Wrapper ICFEP/ABAQUS with UQLink:
modelname = 'NonLinear_Truss';
filename = 'NonLinear_Truss';
switch AnParam.FE_package
case 'ABAQUS'
myUQLinkModel = ABAQUS_wrapper(modelname, filename);
case 'ICFEP'
myUQLinkModel = ICFEP_wrapper(modelname, filename);
end
Part 2: EoD + Bayesian inference
2.2 Recursive inference
for jj = 1:AnParam.N_Obs_stage
% enrich EoD manually and extract enriched FE restults
if jj==1 %First stage EoD run
[~] = EnrichedFERun(myUQLinkModel,Prior_EoD, AnParam,filename,directory,jj);
else
[~] = EnrichedFERun(myUQLinkModel,Prior_EoD, AnParam,filename,directory,jj, Posterior);
end
% Merge different EoD results to EDarray/FEcell_output if possible
% Update the number of FE runs if possible
[coord, EDarray, FEcell_output,AnParam,FE_stageTime] = merge_EoD(AnParam);
%%Note: (1) FE_stageTime and observation_stageTime may be not consistant,
% (2) FE coordination and observation coordination may be not consistent
%%interpolation to make FE result to fit with observation
%%(FE_stageTime/FE coordination must contain observation_stageTime/observationcoord)
FEcell_output_interp = interp_coord_stageTime(observationcoord,coord,FEcell_output, ...
AnParam,observation_stageTime,FE_stageTime);
%% Trim coordination and FE results if possible (trim outputput 1)
[coord,FEcell_output_interp,observationsarray] = Trim(jj,AnParam,FEcell_output_interp, ...
coord,1,observationsarray);
%% Construct PCA-PCE for a given stage {jj-th stage}
[myPCE, PCA] = PCEPCA(AnParam, EDarray, FEcell_output_interp, jj);
%% Make PriorOpts/myPriorDist structure for parameter/observation discrepancy priors
%% Initial stage prior should be defined manually
if jj ==1
%intial srage prior should be defined manually
PriorOpts = []; Posterior{jj} = []; %initial empty cells
[PriorOpts, myPriorDist] = Dynamic_PriorOpts(Prior_forward,jj,AnParam, ...
observationsarray,Prior_hyperdiscrep,Posterior);
else
%rest stages should be inferred automatically
PriorOpts = []; %clear PriorOpts from last stage
[PriorOpts, myPriorDist] = Dynamic_PriorOpts(Prior_forward,jj,AnParam, ...
observationsarray,Prior_hyperdiscrep,Posterior);
end
%%create empty folderes to restore results at initial stage
[~] = emptyFolder(jj)
%%%prior predictive
[~] = predictive_prior(AnParam,EDarray,FEcell_output_interp,coord,observationsarray, ...
gtarray,observation_stageTime,directory,myPriorDist,jj);
%% === Observations ====
%%UQLAB doesnot support NAN in the myData.y
%%Repalce NAN with numerical value eplision=1e-20
epsilon = 1e-20;
for kk = 1:AnParam.N_outputfields
if ~isnan(observationsarray{kk}(:,jj))
obsconcat{kk} = observationsarray{kk}(:,jj)';
else
obsconcat{kk} = epsilon*ones(size(observationsarray{kk}(:,jj)))';
end
end
myData.y = obsconcat;
myData.Name = 'Multiple observations data sources';
% %% === Likelihood switch ====
myLogLikeli = likeliSwitch(jj,observationsarray,PCA,myPCE,AnParam,Prior_hyperdiscrep,coord,directory);
%Add boundary to forward prior (since UQlab not supporting trunction during coupula/marginal input
% inference) we need put additional boundary afterwards
myPriorDist = boundary(PriorOpts,Prior_EoD,AnParam,myPriorDist);
%% Solver options: AIES MCMC
Solver.Type = 'MCMC';
%Solver.MCMC.Visualize.Parameters = 1:AnParam.N_parameters;
%Solver.MCMC.Visualize.Interval = 20;
Solver.MCMC.Sampler = 'AIES';
Solver.MCMC.Steps = 300;
Solver.MCMC.NChains = 100;
BayesOpts.Data = myData;
BayesOpts.Type = 'inversion';
BayesOpts.LogLikelihood = myLogLikeli;
BayesOpts.Solver = Solver;
BayesOpts.Prior = myPriorDist;
BayesAnalysis = uq_createAnalysis(BayesOpts);
%% response surface for RSM
[~] = RMSE_plotting(BayesAnalysis,AnParam,jj,myPCE, ...
PCA,EDarray,observationsarray,gtarray,6000);
%% sobol indices for plotting
[~] = Sobol_indices_plotting(AnParam, FEcell_output_interp,EDarray,jj,observation_stageTime);
%% Postprocess (70% burn in; extract samples from posterior )
Posterior = BI_postprocess(jj,BayesAnalysis,AnParam,Posterior,gtarray,myPriorDist,Prior_forward);
%% save/display MAP for each stage with stages going
[~] = MAP(AnParam,Posterior, directory);
%% plot paramter three sigma range with stages going
[~] = threeSigma_plotting(AnParam,Posterior, gtarray,jj,Prior_forward);
%% Predictive posterior
[~] = predictive_posterior(Posterior,AnParam,jj,EDarray,FEcell_output_interp, ...
coord,observationsarray, gtarray,observation_stageTime,directory);
%% save LL evaluation values at MAP in a csv file if required
[~] = LL_MAP_test(BayesAnalysis,AnParam, jj)
end
